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Abstract 

We have developed a new, very efficient numerical scheme to solve the CR diffusion 
convection equation that can be applied to the study of the nonlinear time evolution 
of CR modified shocks for arbitrary spatial diffusion properties. The efficiency of 
the scheme derives from its use of coarse-grained finite momentum volumes. This 
approach has enabled us, using ~ 10 — 20 momentum bins spanning nine orders 
of magnitude in momentum, to carry out simulations that agree well with results 
from simulations of modified shocks carried out with our conventional finite differ- 
ence scheme requiring more than an order of magnitude more momentum points. 
The coarse-grained, CGMV scheme reduces execution times by a factor approxi- 
mately half the ratio of momentum bins used in the two methods. Depending on 
the momentum dependence of the diffusion, additional economies in required spatial 
and time resolution can be utilized in the CGMV scheme, as well. These allow a 
computational speed-up of at least an order of magnitude in some cases. 
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1 Introduction 



Collisionless shocks are widely thought to be effective accelerators of ener- 
getic, nonthermal particles (hereafter Cosmic- Rays or CRs). Those particles 
play central roles in many astrophysical problems. The physical basis of the 
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responsible diffusive sfiock acceleration (DSA) process is now well established 
through in-situ measurements of heliospheric shocks [10,11] and through ana- 
lytic and numerical calculations [3,6,7,25,31]. While test particle DSA model 
treatments are relatively well developed; e.g., [7,32], it has long been recog- 
nized that DSA is an integral part of coUisionless shock physics and that there 
arc substantial and highly nonlinear backrcactions from the CRs to the bulk 
flows and to the MHD wave turbulence mediating the CR diffusive transport 
(see, for example, [31] and references therein). Most critically, the CRs can cap- 
ture a large fraction of the kinetic energy dissipated across such transitions. As 
they diffuse upstream the CRs form a pressure gradient that decelerates and 
compresses the entering flow inside a broad shock precursor. That, in turn, 
can lead to greatly altered full shock jump conditions, especially if the most 
energetic CRs, which can have very large scattering lengths, escape the system 
and carry significant energy with them. Also in response to the momentum 
dependent scattering lengths and flow speed variations through the shock pre- 
cursor the CR momentum distribution will take on different forms than in a 
simple discontinuity. Effective analytic (e.g., [5,6,28]) and numerical (e.g., [9]) 
methods have been developed that allow one to compute steady-state modi- 
fled shock properties given an assumed diffusion behavior. On the other hand, 
as the CR particle population evolves in time during the formation of such a 
shock the shock dynamics and the CR-scattering wave turbulence evolve as 
well. For dynamically evolving phenomena, such as supernova remnants, the 
time scale for shock modiflcation can be comparable to the dynamical time 
scales of the problem. 

The above factors make it essential to be able to include both nonlinear and 
time dependent effects in studies of DSA. Generally, numerical simulations 
are called for. Full plasma simulations offer the most complete time depen- 
dent treatments of the associated shock microphysics [14,35], but are far too 
expensive to follow the shock evolution over the time, length and energy scales 
needed to model astrophysical CR acceleration. The most powerful alterna- 
tive approach utilizes continuum methods, with a kinetic equation for each 
CR component combined with suitably modifled compressible fluid dynamical 
equations for the bulk plasma (see §2 below). By extending that equation set 
to include relevant wave action equations for the wave turbulence that me- 
diates CR transport, a self-consistent, closed system of equations is possible 
(e.g., [1,2,17,18])). Continuum DSA simulations of the kind just described are 
still quite challenging and expensive even with only one spatial dimension. 
The numerical difficulty derives especially from the very large range of CR 
momenta that must be followed, which usually extends to hundreds of GeV/c 
or beyond on the upper end and down to values close to those of the bulk 
thermal population, with nonrelativistic momenta. The latter are needed in 
part to account for "injection" of CRs due to incomplete thermalization that 
is characteristic of coUisionless shocks. 
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One computational constraint comes from the fact that CR resonant scattering 
lengths from MHD turbulence, \{p), are generally expected to be increasing 
functions of particle rigidity, pc/Ze. The characteristic length coupling the 
CRs of a given momentum, p, to the bulk flow and defining the width of the 
modified shock precursor is the so-called diffusion length, Xd{p) — ^Xv/us, 
where v is the CR particle speed, and Us is the bulk flow speed into the 
shock. One must spatially resolve the modified shock transition for the entire 
range of x^ip) in order to capture the physics of the shock formation and 
the spatial diffusion of the CRs, in particular. The relevant Xd{p) typically 
spans several orders of magnitude, beginning close to the dissipation length 
of the thermal plasma, which defines the thickness of the classical, "viscous" 
gas shock, also called the "subshock" in modified structure. This resolution 
requirement generally leads to very fine spatial grids in comparison to the 
"outer scale" of the problem, which must exceed the largest Xd{p)- 

Two approaches have been applied successfully so far to manage this constraint 
in DSA simulations. Berezhko and collaborators [3] developed a method that 
normalizes the spatial variable by Xd{p) at each momentum value of interest 
during solution of the CR kinetic equation. This approach establishes an spa- 
tial grid that varies economically in tune with Xd{p)- Derived CR distribution 
properties at different momenta can be combined to estimate feed-back on 
the bulk fiow at appropriate scales. The method was designed for use with 
CR diffusion properties known a priori. It is not readily applied to CR diffu- 
sion behaviors incorporating arbitrary, nonlinear feedback between the waves 
and the CRs. As an alternative that can accommodate those latter diffu- 
sion properties, Kang et al.[24] have implemented diffusive CR transport into 
a multi- level adaptive mesh refinement (AMR) environment. The benefit of 
AMR in this context comes from the feature that the highest resolutions arc 
only necessary very close to the subshock, which can still be treated as a dis- 
continuity satisfying standard Rankine-Hugoniot relations. By efficient use of 
spatial gridding both of these computational strategies can greatly reduce the 
cost of time dependent DSA simulations. 

On the other hand, the above methods do not directly address the princi- 
pal computational cost in such simulations, so they remain much more costly 
compared to purely hydrodynamic or MHD simulations. This is because the 
dependence of / on CR momentum, p, adds a physical dimension to the prob- 
lem. In practice, the spatial evolution of the kinetic equation for each CR 
constituent must be updated over the entire spatial grid at multiple momen- 
tum values; say, Np. The value of Np is usually large, since the spanned range 
of CR momentum is typically several orders of magnitude. Physically, CRs 
propagate in momentum space during DSA in response to adiabatic compres- 
sion in the bulk fiow, sometimes by momentum diffusion (see, for example, 
equation 1 below), or because of various irreversible energy loss mechanisms, 
such as Coulomb or radiative losses. The associated evolution rates for / de- 
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pend on the process, but generally depend on df/dlnp. The conventional 
approach to evolving / approximates df/dlnp through low order finite difi^er- 
ences in Alnp (e.g., [3,12]). Experience has shown that converged solutions 
of f{p) using such methods require Alnp < 0.1 [3,21]. In that case, for exam- 
ple, a mere five decades of momentum coverage requires more than 100 grid 
points in In p. Since spatial update of the kinetic equation at each momentum 
grid point requires computational effort comparable to that for any of the ac- 
companying hydrodynamical equations (e.g., the mass continuity equation), 
CR transport then dominates the computational effort by a very large factor, 
commonly exceeding an order of magnitude. 

An attractive alternative approach to evolving the kinetic equation replaces 
/ by its integral moments over a discrete set of finite momentum volumes, 
in which case df/dlnp is replaced by / evaluated at the boundaries of those 
volumes. The method we outline here follows that strategy. Because /(Inp) 
is relatively smooth, simple subvolume models can effectively be applied over 
moderately large momentum volumes. We have found this method to give 
accurate solutions to the evolution of / with an order of magnitude fewer 
momentum bins than needed in our previous finite difference calculations. 
The computational effort to evolve the CR population is thereby reduced 
to a level comparable to that for the hydrodynamics. In recognition of its 
distinctive features we refer to the method as "Coarse-Grained Momentum 
finite Volume" or "CCMV". It extends related ideas introduced in [19], [20] 
and [33] for test particle CR transport. Those previous presentations, while 
satisfactorily following CR transport in many large-scale, smooth flows, did not 
include spatial or momentum diffusion, so could not explicitly follow evolution 
of / during DSA. Instead, analytic, test-particle solutions for f{p) were applied 
at shock jumps. Here we extend the CCMV method so that it can be applied 
to the treatment of fully nonlinear CR modified shocks. 

We outline the basic CGMV method and its implementation in Eulerian hy- 
drodynamics codes in §2. Several tests are discussed in §3, and our conclusions 
are presented in §4. 



2 The Method 



2.1 Basic Equation 

The standard diffusion-convection form of the kinetic equation describing the 
evolution of the isotropic CR distribution function, f{x,p,t), can be written 
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in one spatial dimension as (e.g., [34]). 

— + u—-- (—] ^ + — (k—] + — — (pD^] + S 
dt dx 3 \dx 1 dy dx \ dx J p^ dy \ dy 



where u is the bulk flow speed, y — Inp, k is the spatial diffusion coefficient, 
D is the momentum diffusion coefficient, and 5* is a representative source 
term. We henceforth express particle momentum in units of mc, where m is 
the particle mass. The first RHS term in equation (1) represents "momentum 
advcction'' in response to adiabatic compression or expansion. For simplic- 
ity of presentation equation 1 neglects for now propagation of the scattering 
turbulence with respect to the bulk plasma, which can be a significant infiu- 
ence when the sonic and Alfvenic Mach numbers of the fiow are comparable. 
Although it is numerically straightforward to include this effect, the details 
are somewhat complex, so we defer that to a follow-up work focussed on CR 
transport in MHD shocks. 

Full solution of the problem at hand requires simultaneous evolution of the 
hydrodynamical fiow, as well as the diffusion coefficients, k. and D. Again 
postponing full MHD, the added equations to be solved are the standard 
gasdynamic equations with CR pressure included. Expressed in conservative, 
Eulerian formulation for one dimensional plane-parallel geometry, they are 

d{pu) djfm' + P, + Pe) _ . 

~dr + " ^'^^ 

djpeg) , djpegu + PgU + Pcu) _ dPc 

~dr + - ""^ " 



where Pg and Pc are the isotropic gas and the CR pressure, respectively, Cg = 
Pg/ p(7g — 1) + M^/2 is the total energy density of the gas per unit mass and the 
rest of the variables have their usual meanings. The injection energy loss term, 
L{x,t), accounts for the energy of the suprathermal particles transferred at 
low energy to the CRs. As usual, CR inertia is neglected in such computations, 
since the mass fraction of the CRs is generally tiny. We note for completeness 
that Pc can be computed from / using the expression 

Pmin 



In the simulations described below we set the particle mass, m = 1, for con- 
venience . 
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2.2 The CGVM Method 



As mentioned in §1, the momentum advection and diffusion terms in equation 
1 typically require Ay < 0.1 when using low order finite difference methods 
in the momentum coordinate [3,21]. The resulting large number of grid points 
in y makes finding the solution of equation 1 the dominate effort in simula- 
tions of DSA. On the other hand, previous studies of DSA as well as direct 
observations of CRs in different environments have shown that f{p) is com- 
monly well described by the form / oc p"'^^^^ where q{p) ~ is a slowly 
varying function of y. Thus, we may expect a piecewise powerlaw form to 
provide an efficient and accurate, two-parameter subgrid model for f{y). Two 
moments of f{p) are sufficient to recover the subgrid model parameters. We 
find it convenient to use 



Pi+1 Vi+l 



= / P'f(p)dp^ I exp(3y)/(y)dy, (6) 



and 



Vi 



Pi+1 Ui+i 



j P^f{p)dp^ j exp{4:y) f{y)dy. (7) 



The first of these moments, nj, is proportional to the spatial number density 
of CRs in the momentum bin Ay, = while for relativistic CRs, gi is 

proportional to the energy density or pressure contribution of CRs in the bin. 
Then, for example. 



Qi 



(8) 



where fi = f{pi) = and di = Pi+i/pi with obvious extension to 

Qi. Either of these moments, plus their ratio, gi/piUi, can be used in straight- 
forward fashion (e.g., iteration) to find both fi and the intrabin index, q^. 

To evolve rii and gi we need the associated moments of equation 1 over the 
finite momentum volume bounded by Ay^. The result for rii is 



dui drii du d ( dni\ 



where = {pi + q{Pi)D{pi)/pi)}p'^f{pi) is a fiux in momentum space, with 
p — — p|f|, and where K^i and Sni are averaged over the momentum interval. 
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according to 



and 



Pi+i 



Sm = j p'^S{p)dp. 



Pi 



Extension of the momentum flux term F„. to include other processes such 
as radiative or Coulomb losses is obvious (e.g., [19]). In practice these fluxes 
should be upwinded according to the signs of p and q. Evaluation of -F„- and 
q{pi) at the boundaries of the included momentum range requires application 
of suitable boundary conditions, of course. We usually have set un^+i — 0. 
In most cases we pick a sufficiently large maximum momentum, PNp+i, that 
this condition is important only late in the simulation, if at all. Appropriate 
conditions at the lowest momenta can be more involved, depending on how 
one intends to connect the CR particle distribution to the thermal particle 
distribution, as in the injection models discussed below. 

The g'j-associated moment of equation 1 is 



where Fg. ^ PiF^^^, 



gi{D/p'^)- = Jp!^"^^ pDfdp, and 5'^. is given by an analogous expression to 
equation 11. 

We note that momentum binning in the CGMV scheme is quite flexible, so 
that it can be easily adapted to either uniform or nonuniform momentum 
bin sizes, or to a momentum range that evolves during the simulation. We 
have successfully implemented both nonuniform and evolving momentum bin 
structures, although, for brevity we do not illustrate them here. 

To update the distribution function we simultaneously integrate equations (9) 
and (12) over the timestep At''. Our implementation applies these methods 
in Eulerian hydrodynamical codes, so we split the update into two parts. 
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Spatial advection is carried out by a second order, explicit van Leer scheme. 
The remaining terms are evolved using a second-order, semi-implicit Crank 
Nicholson scheme. For example, given values at timestep t'^, and assuming 
for illustration a uniform spatial grid, the implicit part of the solution is given 
by the tridiagonal system 

Atn^tl + + = Cf, (14) 



where 



^° = l + 5(i^,+l/2 + i^,-l/2), 



At" 



Pi 



with 5 = (l/2)AiV(Ax)2, and i^^+i/s = (l/2)(i^i+i + K^). 

The coefficients in equation 14 are obtained with the aid of the solutions to 
equations 2-4, which are updated prior to solution of equations 9 and 12. We 
note again that similar methods can be applied to follow the evolution of the 
wave turbulence that resonantly scatters CRs and that defines the spatial and 
momentum diffusion coefficients. In that case one begins with the wave action 
equation for the appropriate waves rather than the particle kinetic equation 
(e.g., [2]). 

The solution of equation (14) for either rii or is quite analogous to our 
previous FD methods. Thus, since the CGMV method evolves two quantities 
rather than one, the relative effort required for a given Np is roughly twice 
in the CGMV scheme that required in the FD scheme. Our tests confirm this 
expected scaling. On the other hand, since Np can be dramatically reduced in 
a CGMV simulation the method can still be more efficient by a large factor. 



2.3 A Flux Fraction Injection model 



The most common source term represented by S in equation (1) is injection at 
the shock of low energy CRs from the thermal plasma. There is presently no 
generally accepted theory for that process. However, we have implemented two 
commonly used models successfully into the CGMV scheme. For completeness 
we outline those here. 
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The simplest and one of the most frequently applied injection models as- 
sumes that a small, fixed fraction of the thermal particle flux through the gas 
subshock, Cinj, is injected at a momentum pi^j = Q^Cs2, where a is a constant 
greater than unity and Cg^ is the plasma sound speed immediately downstream 
of the subshock (e.g., [21]). This gives S{p) — €inj{pi/ nmp)usw{x — Xs)S{p — 
Pinj), where pi is the plasma mass density just upstream of the subshock, firup 
is the plasma mean particle mass, Us is the subshock speed with respect to 
the plasma immediately upstream and w is a normalized weight function that 
allows the injection to be distributed across the numerical shock structure. 
In this case = ^einj{pi/ iJ,mp)usw{x - x^), while Sg^ = Pmj-S'„. in the mo- 
mentum bin with pi < pinj < Pi+i- Both Sm and Sg^ are zero, otherwise. The 
energy extracted from the thermal plasma is simply L = ^einjw{x)a'^c^^piUs- 
For convenience we call this injection model the "flux fraction" or "FF" model. 

2.4 A Thermal Leakage Injection model 

A more sophisticated approach to injection physics includes models of the 
physical processes moderating particle orbits in the post shock flow region in 
order to estimate the probability that particles of a given speed will be able to 
escape back upstream, across the subshock. In such "thermal leakage" (TL) 
models for CR injection at shocks, most of the downstream thermal protons 
are locally confined by nonlinear MHD waves and only particles well into 
the tail of the postshock Maxwellian distribution can leak upstream across 
the subshock. In particular, "leaking" particles not only must have velocities 
large enough to swim against the downstream flow in order to return across the 
shock, they must also avoid being scattered during that passage by the MHD 
waves that mediate the plasma subshock. To model TL injection we utilize a 
"transparency function", Tesc, expressing the probability that supra-thermal 
particles at a given velocity can leak upstream from behind the subshock (see 
[15] for details). In particular we set 

(1 - -) 

Tesciv, U2) =H[d-{l + cb)] exp {- [d-{l + 63)]-'} , (15) 

where U2 is the postshock flow speed in the subshock frame, H is the Heaviside 
step function, and the particle velocity is normalized to v = veB/u2. The 
parameter, = Bq/B±, measures the ratio of the amplitude of the postshock 
MHD wave turbulence B± to the general magnetic field aligned with the shock 
normal, Bq. Both hybrid simulations and theory suggest that 0.25 < 
0.35 [30], so that the model is well constrained. With this Tgsc the shock is 
completely "opaque" to particles with momenta less than pi, i.e., Tgsc = 
for p < pi, where pi = mpU2{l + eB)/eB- So pi(t) is the lowest momentum 
of the first momentum bin in the TL model and changes in time with the 
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postshock flow speed. For sb ~ 0.3, pi ~ 4 — ^mipiia- The shock becomes 
virtually transparent to particles with momenta two to three times greater 
than pi. For strong, unmodified shocks pi in the TL model and Pinj in the FF 
model are similar when eb ~ 0.3 and a ~ 2. Under those circumstances the 
initial injection rates will be roughly similar, although differences in model 
physics lead to different behaviors as such shocks become modified (see, e.g., 
[4], [15], [25]). 

The TL model is implemented in the CGMV scheme by the following numer- 
ical approach. After solution of equations (9) and (12) the net changes in rii 
and Qi are corrected (reduced) in the upstream region by application of the 
transparency function as follows: 

n1'-'^v!l+ j TUp){ft^' - ft)p'dp (16) 

Pi 



and 

Pi+1 

9t^'-9- + J resc{p){fT' - f')p'dp (17) 
Pi 



where f^^^, found using equation (8), is the CR distribution updated with 
equations (9) and (12). The energy loss rate of the bulk plasma to injection 
into the i-th CR momentum bin can be approximated by 

Hx,t) ^ -^^{^ 7 ^P"(^^^-^)f^dP (18) 

Pi 



(see [25]). 

With the piece- wise power-law subgrid model {fi{p) = fi{p/Pi)~'^') the inte- 
grals in equations (16)-(18) can be written: 



7 rap)Mp)p'dp = 7 ^-cP^^-^^^rfi/, (19) 

Pi 



riMi - 3)p! 
(1 - d^' 

(1 - 

Pi ^ ^ ' Pi 



j nsc{p)Mp)p'dp = ^^'j/^' / rescP^'-'-^dy, (20) 
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and 



In practice, this leakage step is significant only for the lowest few momentum 
bins, so that this correction need not be applied to all bins. 



3 Discussion 

In order to test the performance of the new CGMV scheme we have installed it 
into two distinct one-dimensional Eulerian hydrodynamic (HD) codes that we 
have previously applied to studies of CR-modified shocks using conventional 
finite difference (FD) methods to solve the diffusion convection equation. In 
this section we briefiy discuss the results and compare the CGMV and FD 
behaviors. Both of the host HD codes are constructed from high order, con- 
servative Riemann solver-based schemes designed to capture shocks sharply. 
First we describe results from the CGMV scheme installed in a second order 
"Total Variation Diminished" (TVD) HD code based on the finite difference 
scheme of Hartcn [16]. This is the HD version of the MHD-CR code used by 
us in a previous study of CR modified shocks [23]. Gas subshocks in the TVD 
scheme typically spread over 2-3 numerical zones. An outline of the code me- 
chanics and the FD CR scheme can be found there, in [13,21] and references 
cited in those papers. The FD solver employed a Crank- Nicholson routine orig- 
inally introduced in [12] for evolving p'^f{p) that is similar to equation (14). 
For the TVD tests we applied the FF injection model. 

In addition we present CGMV tests carried out with our Cosmic Ray Adap- 
tive refinement SHock (CRASH) code. CRASH is based on the high order 
Godunov-like shock tracking algorithm of LeVeque [27]. The hydrodynam- 
ics routine in that code employs a nonlinear Riemann solver to follow shock 
discontinuities within the zones of an initially uniform grid. Thus, gas sub- 
shocks in CR-modified shocks remain discontinuous throughout a simulation, 
allowing CR transport to be modeled down almost to the scale of the physi- 
cal shock thickness. CRASH also employs adaptive mesh refinement (AMR) 
around shocks in order to reduce the computational effort on the spatial grid. 
Refinement is centered on the subshock and each level spans 100 zones with a 
resolution twice as fine as the level above it. The number of refinement levels 
depends on what is required to capture diffusion of the lowest energy CRs. 
The standard, previously documented version of CRASH uses the same FD 
methods as the TVD code to solve the diffusion convection equation for CR 
transport. It is described in [24] and [25]. For the CRASH tests discussed in 



11 



this paper we employed the TL injection model. 



3.1 TVD-CR Tests 



We first examine some results obtained using the TVD-CR code with both 
FD and CGMV schemes used to model the evolution of a strong CR-modified 
shock. Fig. 1 and Fig. 2 illustrate the evolution of shocks formed by the reflec- 
tion of a Mach 30 flow (adiabatic index, 7 = 5/3) off the left grid boundary. 
The resulting piston-driven shock initially has a Mach number, Mg ~ 40. The 
density and sound speed of plasma entering from the right boundary were 
set to p = 1 and Cg = respectively, so that the inflow speed was unity. 
The time unit for the calculations is also set by these scalings. In order to 
relate hydro dynamical variables to CR momenta it is necessary to fix the 
unit flow speed (the inflow speed in all the simulations discussed in this pa- 
per) with respect to the speed of light; i.e., c = l/f3. We set f3 = 10~^ in the 
TVD-CR simulations. Time steps were flxed by a standard Courant condition. 
At = 0.8Ax/max{u±Cs). 

For the simulation illustrated in Fig. 1 evolution of the CR distribution is 
followed over the momentum range [pmin,Pmax] = [2 x 10"'^, 1.6 x lO''] {y„iax — 
Umin ~ 16). The simulation represented in Fig. 2 included the momentum 

range \pmin,Pmax\ = [2 X 10""^, 2.4 X 10^] {ijmax - Vmin ~ 21). 

The CR diffusion coefficient, n is spatially uniform and set to K{p) = Q.lp^'^^. 
In a quest for a reasonably generalized behavior that required minimal compu- 
tational effort, this choice was motivated by results from Malkov, who found 
self-similar analytic steady-state solutions for strong CR-modified shocks that 
apply to all powerlaw forms of /«(p), so long as the powerlaw index is steeper 
that 0.5 [29]. Thus, our k choice leads to fairly general shock behaviors in a 
way that minimizes the width of the precursor. That width, which determines 
the minimum space that must be simulated ahead of the subshock, is set by 
Xd{Pcutoff), where Pcutoff represents the maximum momentum contained. The 
TVD-CR simulations utilize a uniform, fixed grid, so, for example, the Bohm 
diffusion form modeled in the CRASH simulations below would lead to ex- 
cessive costs for the TVD-CR FD test simulations presented in this section. 
The spatial resolution required for the calculations is set by Xd{pmin), since 
accurate solutions of the diffusion-convection equation require good structural 
information in the diffusive shock precursor upstream of the subshock. Pre- 
vious convergence tests have shown that Ax < O.lxd is desirable (e.g., [21]). 
For the problems illustrated in Fig. 1 and Fig. 2 these considerations led us 
to set Ax = 3.8 x 10"^ for both the FD and CGMV simulations. By varymg 
this resolution, we verified that the shock evolution is reasonably converged 
with respect to Ax. 
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The simulation followed in Fig. 1 assumes a pre-existing CR population, 
f{p) 0*^ corresponding to an upstream CR pressure. = Pg. No fresh 

injection is included at the shock; i.e., einj = 0. This test then provides a sim- 
ple and direct comparison between the CGMV and FD schemes for solving 
the diffusion-convection equation, since it omits any complications related to 
the injection model. 

This simulation pair evolves the shocks until t = 30, which is sufficient to 
accelerate CRs to ultrarelativistic momenta. The spatial grid spans the interval 
[0,16], which is sufficient to contain the leading edge of the CR precursor to 
the end of the simulation. 

The CR- modified shock spatial structures and the CR momentum distribu- 
tions at the subshock are shown in Fig. 1 at times t — 2, 10, 20, and 30. Before 
comparing the solutions obtained with the two different methods, it is useful 
to summarize briefly the physics captured during the shock evolution. All the 
behaviors described here have been reported previously by multiple authors. 
The figure shows the well-known property of strongly modified shocks that 
the CRs extract most of the energy flux into the structure. That leads to a 
substantial drop in the postshock gas pressure, Pg, and a large increase in the 
postshock density, p. Together those indicate a strongly reduced postshock 
gas temperature. The decreased temperature is evident in the p^/ plot at the 
subshock, which is dominated at low momenta by the Maxwellian distribution 
of the bulk plasma. As CRs diffuse upstream against the inflowing gas they 
compress the flow within the precursor, preheating the gas (adiabatically in 
these simulations). Initially, while the CR pressure is relatively small compared 
to the incoming momentum flux, the gas subshock remains strong enough to 
produce a full four times density compression on top of the prccomprcssion. 
However, the subshock weakens once Pc ~ pu^, reducing the subshock com- 
pression in this case to a factor 2.6, corresponding to a subshock Mach 
number near 2.3. That evolution explains the well-known, transitory "density 
spike" in the shock structures seen after t = 2. We note that since energy 
extracted from the flow by CRs becomes increasingly spread upstream and 
downstream due to CR diffusion, the total compression in such an evolving 
modified shock would always exceed the factor of seven one would predict 
for a strong, fully relativistic gas shock. For this simulation no significant CR 
energy escapes the spatial grid through upstream diffusion. However, at late 
times (t > 20) the partial pressure due to CRs just below Pmax is sufficient 
that escape across the upper momentum boundary is significant. This con- 
tributes to the slow decrease in Pc behind the shock and the increasing total 
compression through the shock transition that is visible in Fig. 1. 

In the early evolutionary stages of this flow, while shock modiflcation is mod- 
est, the CR momentum distribution resembles the powerlaw form, p'^f{p) — 
const, predicted by test particle theory for a strong shock (e.g., [7]). With the 
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spatial diffusion coefficient used in tfiis simulation the high momentum cutoff 
to the distril^ution increases with time approximately as Pcutoff (see, e.g., 
[26]). As shock modification intensifies, most of the fiow compression shifts 
from the subshock to the precursor. Then DSA of high momentum CRs oc- 
curs predominantly within the precursor rather than near the less important 
subshock. Consequently the CR distribution develops the familiar upwards- 
concave form resulting from the momentum dependent CR diffusion length. 
CRs of higher momentum experience a greater velocity jump within the pre- 
cursor, so gain more energy each time they are reflected within the shock 
structure. That flattens the distribution, f{p) at momenta below Pcutoff- The 
result is a bump in the distribution of p'^f{p)- On the other hand, CRs with 
momenta only a little above the injection range remain trapped close to the 
subshock. Their distribution closely approaches the steady state, powerlaw, 
test particle form appropriate to the weakened subshock. That feature ex- 
tends upwards in momentum as the bump near p cutoff moves upwards. 

Looking finally to compare the two methods used to evolve the shock evolu- 
tion displayed in Fig. 1 we see results from the FD scheme with Ay = 0.11 
and the CGMV scheme with Ay — 1.0. The agreement is generally very good. 
All the dynamical quantities, including shock jumps and the CR momentum 
distributions show excellent agreement. The curves representing the FD and 
CGMV distributions of p, Pg and Pc and virtually indistinguishable in the 
plots. Most notably, all features formed in the FD evolution of the CR mo- 
mentum distribution are faithfully reproduced by the much coarser CGMV 
distribution. 

Given the excellent comparison in this strongly modified fiow it is satisfying to 
note that the execution time required for the CGMV solution was a little less 
than 20% of that for the FD solution, demonstrating the significantly higher 
efficiency of the former method. The speed-up observed in our implementations 
of the two methods is roughly in accordance with what we would predict for 
a given reduction factor in the number of momentum values used, since the 
CGMV method requires one to evolve two distributions rii and Qi for each 
momentum bin. We address convergence with respect to momentum resolution 
in our discussion of a second shock. 

Fig. 2 shows a Mach 30 fiow similar to that in Fig. 1. In this case FF injection 
is included with commonly assumed values, einj = 10^^ and a = 2.0 (sec §2.3), 
the upper momentum bound is increased to Pmax = 2.4 x 10^ and the spatial 
grid extends farther from the piston to x^ax — 25. A negligible pre-existing 
CR population is included to avoid numerical issues coming from the fact that 
our CGMV scheme requires computation of the ratio gi/piUi over the entire 
grid. 

The simulations evolved the shock until t — 70, which leads to Pcutoff ^ 10^. 
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The spatial grid, spanning the interval [0,25], is sufficient to contain the CR 
leading edge of the precursor almost to the end of the simulations. However, 
after t ~ 50 some CR energy escapes through the right boundary, due to 
diffusion upstream, mimicking the behavior of a "free escape" boundary (FEB) 
(e.g., [22] and references therein). Just as for the shock simulated in Fig.l, this 
energy loss amounts to a cooling process, so that the total shock compression 
increases with time as the simulation ends. 

Again the agreement between FD and CGMV simulations is very good. Very 
early in the evolution of the shock, when the CRs are dominated by freshly 
injected nonrelativistic particles, the shock evolution is slightly faster in the 
CGMV scheme. That influence becomes insigniflcant later on, so that the 
modified shock structures found by the two schemes are almost identical as is 
the distribution p^f{p) at the shock. There is a small residual effect that the 
position of the subshock in the CGMV simulation lags slightly behind that 
of the FD simulation, and that Pc is slightly higher in the CGMV simulation 
near the piston, where the shock first formed. 

Since the efficiency of the CGMV method comes from its ability to cover 
the momentum range coarsely, it is important to evaluate how broad the 
momentum bins can be and still faithfully model the evolution of the shock. 
Fig. 3 illustrates convergence of the CGMV scheme with respect to momentum 
bin size. Ay, at t = 70 for the flow modeled in Fig. 2. The upper panel plots the 
spatial Pc distributions, while the lower panel shows the particle momentum 
distributions at the gas subshock. For reference the corresponding FD solution 
(Ay = 0.11) is shown by the dotted red curves. Solutions from the CGMV 
scheme are plotted for Ay = 1, Ay = 1.40 and Ay = 2.1. Even the coarsest 
of the CGMV sohitions is in basic agreement with the other CGMV solutions 
and with the FD solution. Fine details in the momentum distribution are 
naturally obscured as the CGMV bin size increases. The largest bins with 
Ay = 2.3 span a decade in CR momentum (pi+i/pi — 10.2), but still capture 
the basic dynamical properties of the CRs correctly. 

However, the quality of the CGMV solutions deteriorates for still larger mo- 
mentum bins in these experiments, once the simple subgrid model for the 
momentum distribution becomes inadequate. As illustrated in the lower panel 
of Fig. 3 already momentum bins larger than roughly Ay ~ 1.5 cannot closely 
follow sharper structures in the momentum distibution that develops at the 
ends of the CR distribution. That enhances Pc upstream of the subshock, 
where the flow is both cold and strongly compressed as it approaches the 
subshock. When the errors become excessive, for Ay > 2.3, in this case, Pc- 
induced overcompression in this region can cause the Riemann solver in our 
TVD code to perform poorly or even to fail in high Mach number flows. In 
general the largest allowed momentum bin size, Amax should depend on the 
strongest curvature of the CR momentum distribution function as well as the 
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degree of shock modification. 



3.2 CRASH Tests 



For a third test example we illustrate in Figs. 4 and 5 simulation results using 
the CRASH code with the TL injection model applied to a Mach 10 flow 
reflecting off the left computational boundary. The initial gas shock Mach 
number is approximately 13. As for the previous test, the upstream gas density 
and flow speed are set to unity, with upstream sound speed, Cg — and 
f3 — 5 X 10~^. The time unit is defined accordingly. CR momenta are tracked 
over the range [pi, W^], where pi{t) is the smallest momentum that can leak 
upstream (see equation 15). 

In this case a Bohm-type diffusion model with kb — p'^/Vp^ + 1) is adopted 
and the TL injection parameter, = 0.2 is used. The CRASfl test was 
significantly more computationally demanding than the TVD-CR tests. Note 
first that in the CRASH simulation the value of kb{p = 1) = l/\/2, while 
k{p = 1) = 0.1 in the previous examples shown in Figs. 1-3. Consequently, 
X(i{p — 1) is about seven times greater in the current case, and the nomi- 
nal physical scale of the precursor and its formation timescale are similarly 
lengthened. In addition, the stronger momentum dependence of Bohm dif- 
fusion coefficient means that the precursor width expands more strongly as 
Pcutoff increases. The associated time rate of increase in Pcutoff is, however, 
slower, so that the shock must evolve longer to reach a given Pcutoff- These 
factors substantially increase the size of the physical domain needed to reach 
a given Pcutoff- 

Fig. 4 shows the early evolution of this CR-modified shock for i < 20 as 
computed with both the FD and the CGMV methods. The spatial domain 
for this simulation is [0,20]. The base spatial grid included 10'^ zones, giving 
Axq = 2 X 10^'^. Since it is necessary to resolve structures near the subshock 
on scales of the diffusion length for freshly injected, suprathermal CRs, the 
AMR feature of the CRASH code is utihzed. The FD simulation is carried out 
with 7 refined grid levels; four levels of refined grid are applied in the CGMV 
simulation. 240 momentum points {Ay = 0.058) are used in the FD simulation, 
while the CGMV simulation includes 20 momentum bins {Ay = 0.72). The 
time step for each refinement level, Ati^ , is determined by a standard Courant 
condition, that is, Ati^ — 0.3 Axi^/ msix{u±Cs). Although the Crank-Nicholson 
scheme is stable with an arbitrary time step, the diffusion convection equation 
is solved with the time step smaller than Atocig ~ '^Ay{Axig/us) to maintain 
good accuracy in the momentum space advection {i.e., dy/dt = — ||f)- With 
Ay = 0.058, the required time step is smaller by a factor of three or so than the 
hydrodynamic time step in the FD simulation. Consequently, the FD diffusion 
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convection solver is typically subcycled about 3 times with Atocig for each 
hydrodynamic time step. Because of the much larger Ay. subcycling is not 
necessary in the CGMV simulation. That adds another relative economy to 
the CGMV calculation. 

At the end of this simulation, t — 20, the modified shock structure is ap- 
proaching a dynamical equilibrium in the sense that the postshock values of 
p, Pg and Pc will not change much at later times. Since this shock is weaker 
than the Mach 40 shocks examined earlier modifications are more moderate. 
On the other hand, as expected from the stronger momentum dependence of 
Xd{p), the shock precursor broadens much more quickly in the present case. 
The cutoff in the CR distribution has reached roughly Pcutoff ~ 10 by i = 20. 
Longer term evolution of this shock will be addressed below. 

The agreement between the FD and CGMV solutions shown in Fig. 4 is good, 
although not as close as it was in the examples illustrated in Fig. 1 and Fig. 2. 
The more apparent distinctions between the two solutions in the present case 
come from effective differences in the application of the TL injection model 
with Bohm diffusion in FD and CGMV methods. Recall that the CGMV 
scheme applies the diffusion coefficients averaged across the momentum bins 
(see equations 10, 13). The Bohm diffusion model has a very steep momen- 
tum dependence for nonrelativistic particles; namely, k oc p^. At low momenta 
where injection takes place the averaging increases the effective diffusion coef- 
ficient, and, thus, the leakage flux of suprathermal particles, leading to higher 
injection rate compared to the FD scheme for the same TL model parameters. 
Consequently, the distribution function in the second bin at p2 ~ 1-4 x 10~^ is 
slightly higher in the CGMV scheme, as evident in Figs. 4-5. Note that f{pi) 
is anchored on the tail of Maxwcllian distribution. The CGMV solutions ac- 
cordingly show slightly more efficient CR acceleration than the FD solutions 
at early times. In this test Pc is about 5 % greater in the CGMV simulation 
at t — 20. Since the CGMV scheme can be implemented with nonuniform 
momentum bins, such differences could be reduced by making the momentum 
bins smaller at low momentum in instances where the details relating to the 
injection rate were important. 

We show in Fig. 5 the evolution of this same shock extended to t = 10"^, as com- 
puted with the CGMV method. This simulation is computed on the domain 
[0,800], spanned by a base spatial grid of 2 x 10^ zones, giving Axq = 4 x 10^^. 
We also included 7 refined grid levels at the subshock, giving Ax^ = 3.1 x 10^"^. 
This grid spacing is insufficient for convergence at the injection momentum, 
Pinj ~ 10~^, so that the very early evolution is somewhat slower than in 
the simulations shown in Fig. 4. However, once shock modification becomes 
strong evolution becomes roughly self-similar, as pointed out previously [25]. 
The time asymptotic states do not depend sensitively on the early injection 
history. The self-similar behavior results with Bohm diffusion from a match 
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between the upstream and downstream extensions of the CR population. One 
also sees from the form of the distribution function in Fig. 5 that the post- 
shock gas temperature has stabilized, while the previously-explained concave 
form to the CR distribution is better developed than it was at earher times. 

This simulation illustrates nicely the relative efficiency of the CGMV scheme. 
The equivalent FD simulation would be very much more expensive, because 
this model requires a long execution time and a large spatial domain. With 
Bohm diffusion k oc p for ultrarelativistic CRs, so that the scale of the precur- 
sor, Xd oc Pcutoff- At the same time the peak in the CR momentum distribution 
extends relatively slowly, with Pcutoff oc t. The required spatial grid is, thus, 
40 times longer than for the shorter simulation illustrated in Fig. 4. The sim- 
ulated time interval in the extended simulation was 50 times longer. Together 
those increase the total computational time by a factor 2000. The FD calcula- 
tion with X — [0, 20] to t = 20 took about 2 CPU days on our fastest available 
processor, so the extended simulation would have been unrealistic using the 
FD method. The extended CGMV simulation, however, required only about 
10 times the effort of the shorter FD simulation, clearly demonstrating the 
efficiency of the CGMV scheme. This speed-up is a result of combination of 
several factors: 20 times larger grid spacing, no need for subcychng for the 
diffusion convection solver, and, of course, a smaller number of momentum 
bins. 



4 Conclusions 

Detailed time dependent simulations of nonlinear CR shock evolution are very 
expensive if one allows for inclusion of arbitrary, self-consistent and possibly 
time dependent spatial diffusion, as well as various other momentum depen- 
dent transport processes. The principal computational cost in such calcula- 
tions is typically the CR transport itself, and, in self-consistent calculations, 
the analogous transport of the MHD wave turbulence that mediates CR trans- 
port. Tracking these behaviors requires adding at least one physical dimension 
to the simulations compared to the associated hydrodynamical calculations, 
since the coUisionless media involved are sensitive to the phase space config- 
urations of the particles and waves. Particle kinetic equations (commonly the 
so-called diffusion convection equation) provide a straightforward approach to 
addressing this problem and can be coupled conveniently with hydrodynamical 
equations that track mass and bulk momentum and energy effectively. 

Momentum derivatives of the CR distribution function in the diffusion convec- 
tion equation are most frequently handled by finite differences. Although it is 
simple, that approach requires moderately fine resolution in momentum space. 
That is a primary reason that such calculations are costly. Here we introduce a 
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new scheme to solve the diffusion convection equation based on finite volumes 
in momentum space with a momentum bin spacing as much as an order of 
magnitude larger than that of the usual finite difference scheme. We demon- 
strate that this Coarse Grained Momentum finite Volume ( CGMV) method 
can be used successfully to model the evolution of strong, CR-modified shocks 
at much lower computational cost than the finite difference approach. The 
computation efficiency is greatly increased, not only because the number of 
momentum bins is smaller, but also because the required spatial grid spacing 
is less demanding due to the coarse-grained averaging of the diffusion coeffi- 
cient used in the CGMV method. In addition, larger momentum bin size can 
eliminate the need of subcycling of the diffusion convection solver that can be 
necessary in some instances using finite differences in momentum. Thus, the 
combination of the CGMV scheme with AMR techniques as developed in our 
CRASH code, for example, should allow more detailed modeling of the diffu- 
sive shock acceleration process with a strongly momentum dependent diffusion 
model such as Bohm diffusion, or self-consistent treatments of CR diffusion 
and wave turbulence transport. 
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Fig. 1 . Evolution of a CR modified Mach 40 plane shock using two different schemes 
for solving the diffusion convection equation. The lower-left panel illustrates the 
momentum distribution at the shock. The other panels show the spatial distribution 
of key dynamical variables. Black dotted lines represent solution with a conventional 
finite difference scheme using Ay = 0.11. red lines and 'stars' were obtained using 
the new CGMV scheme with Ay = 1.0. The solutions are almost indistinguishable. 
A pre-existing CR population, f{p) oc corresponding to the upstream CR 

pressure, Pc = Pg is included, without fresh injection at the shock (ej„j = 0.). 
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Fig. 2. Same as Fig. 1 except the flux fraction injection model is adopted with no 
pre-existing CR population. The momentum distribution includes the Maxwellian 
distribution of the thermal plasma based on its temperature. 
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in Fig 2. The different curves represent results computed using the FD scheme and 
three different momentum resolutions with the CGMV scheme. Bottom: The CR 
distribution function at the shock from the same simulations. 
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Fig. 4. Evolution of a CR modified Mach 13 plane shock using the CRASH code 
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Fig. 5. Same as Fig. 4 except that the terminal time and simulation box are ex- 
tended to t = 10'^ and Xmax = 800, respectively. The heavy dashed lines represent 
solution at t = 20 with a conventional finite difference scheme using 240 momentum 
points {Ay = 0.058). The red solid lines and 'X's represent CGMV solutions at 
t = 20,200,600, and 1000 with 20 momentum bins (Ay = 0.72). 
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